import matplotlib.pyplot as plt
    from scipy.io import wavfile
    import numpy as np
    from IPython.display import Audio
    import scipy.signal as sig
    import scipy.fft as fft

%matplotlib widget

def fm_modulation(
        modulating_signal: np.ndarray,
        carrier_freq: float,
        fs: float,
        mod_index: float
) -> np.ndarray:
    """
    Frequency modulation
    Parameters
    ----------
    modulating_signal: numpy.ndarray
        modulating signal
    carrier_freq: float
        carrier frequency
    fs: float
        sampling frequency
    mod_index: float
        modulation index

    Returns
    -------
    numpy.ndarray
        modulated signal
    """
    dt = 1/fs
    time = np.arange(len(modulating_signal)) * dt
    modulated_signal = np.sin(2 * np.pi * carrier_freq * time + mod_index * np.cumsum(modulating_signal)*dt)
    return modulated_signal

def am_demodulation(
        modulated_signal: np.ndarray,
) -> np.ndarray:
    """
    Amplitude demodulation
    Parameters
    ----------
    modulated_signal: numpy.ndarray
        modulated signal
    fs: float
        sampling frequency
    mod_index: float
        modulation index

    Returns
    -------
    numpy.ndarray
        demodulated signal
    """
    # use hilbert transform to get the envelope
    analytic_signal = sig.hilbert(modulated_signal)
    amplitude_envelope = np.abs(analytic_signal)
    # demodulate
    return amplitude_envelope

def fm_demodulation(
        modulated_signal: np.ndarray,
        carrier_freq: float,
        fs: float,
        mod_index: float
) -> np.ndarray:
    """
    Frequency demodulation
    Parameters
    ----------
    modulated_signal: numpy.ndarray
        modulated signal
    carrier_freq: float
        carrier frequency
    fs: float
        sampling frequency

    Returns
    -------
    numpy.ndarray
        demodulated signal
    """
    # differentiate the modulated signal and than use am demodulation to get the instantaneous phase
    dt = 1/fs
    diff_modulated_signal = np.diff(modulated_signal) / dt
    demod = am_demodulation(diff_modulated_signal)
    return demod - np.mean(demod)

Solution

Audio("fm_modulated_fs_250khz_fc_80khz.wav")

    fs, data_mod = wavfile.read("fm_modulated_fs_250khz_fc_80khz.wav")
    fc = 80e3
    Audio(fm_demodulation(data_mod, fc, fs, 10000), rate=fs)

make the EX

    # load the sound file
    fs, data = wavfile.read("piano_250khz.wav")
    data = data / np.max(np.abs(data))
    data = data - np.mean(data)
    Audio(data, rate=fs)

freq_axis = fft.fftshift(fft.fftfreq(len(data), 1/fs))
data_f = fft.fftshift(fft.fft(data))
plt.figure()
plt.plot(freq_axis, 20*np.log10(np.abs(data_f)))

[<matplotlib.lines.Line2D at 0x20d4cea37c0>]

2**15

32768

data_mod.min()

np.int16(-1023)

fc = 80e3
data_mod = fm_modulation(data,fc, fs, 20000)
data_mod = data_mod * (2**10)
data_mod = data_mod.astype(np.int16)
data_mod_f = fft.fftshift(fft.fft(data_mod))
plt.figure()
plt.plot(freq_axis, 20*np.log10(np.abs(data_mod_f)))

[<matplotlib.lines.Line2D at 0x20d4e56b760>]

np.save("fm_modulated_fs_250khz_fc_80khz.npy", data_mod)
wavfile.write("fm_modulated_fs_250khz_fc_80khz.wav", fs, data_mod)

Audio("fm_modulated_fs_250khz_fc_80khz.wav")

Audio(fm_demodulation(data_mod, fc, fs, 10000), rate=fs)

K = int(4*fs)
iq_data = data[:K,0] + 1j*data[:K,1]

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[14], line 2
      1 K = int(4*fs)
----> 2 iq_data = data[:K,0] + 1j*data[:K,1]

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

plt.plot(iq_data.real[:500])
plt.plot(iq_data.imag[:500])
plt.show()

fc = 80e3
time_axis = np.arange(len(iq_data)) / fs
iq_passband = iq_data * np.exp(2j*np.pi*fc*time_axis)
iq_air = np.real(iq_passband)

my_demod = fm_demodulation(iq_air, fc, fs, 1.0)
display(Audio(my_demod, rate=fs))

502654.8245743669

iq_demode = fm_demod(iq_data, df=1.0, fc=0.0)
display(Audio(iq_demode, rate=fs))

iq_data_modulated = np.real(iq_data)

def fm_modulation(
        modulating_signal: np.ndarray,
        carrier_freq: float,
        fs: float,
        mod_index: float
) -> np.ndarray:
    """
    Frequency modulation
    Parameters
    ----------
    modulating_signal: numpy.ndarray
        modulating signal
    carrier_freq: float
        carrier frequency
    fs: float
        sampling frequency
    mod_index: float
        modulation index

    Returns
    -------
    numpy.ndarray
        modulated signal
    """
    dt = 1/fs
    time = np.arange(len(modulating_signal)) * dt
    modulated_signal = np.sin(2 * np.pi * carrier_freq * time + mod_index * np.cumsum(modulating_signal)*dt)
    return modulated_signal

fc = 80e3
modulated_signal = fm_modulation(iq_demode, fc, fs, 1.0)

display(Audio(fm_demod(modulated_signal,1, 2*fc/fs), rate=fs))